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ABSTRACT 


This thesis analyzes the sway, yaw, and roll dynamic stability of neutrally 
buoyant submersibles. Utilizing the hydrodynamic coefficients for a Mark 
IX Swimmer Delivery Vehicle (SDV) as a base-line model, the linearized 
equations of motion for the decoupled steering and roll systems are 
compared to the coupled system. Two different configurations of 
hydrodynamic coefficients are considered along with the effects of varying 
the vertical (Z,) and longitudinal (X,) centers of gravity of the vehicle while 
the longitudinal center of buoyancy (Xp) is held constant. Results 
demonstrate the significant effects on stability of coupling the steering and 
roll equations of motion, and the importance of Z,; and Xg selection in 
minimizing those effects while retaining stability. Perturbation analysis 
results confirm the essential dependence of the linearized coupled 


equations on the separation of X, and Xp. 
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I. INTRODUCTION 


A. GENERAL 


As the missions for submersibles become more complex and 
demanding, the requirement for a highly stable platform becomes 
increasingly important so the operator(s) can concentrate on matters other 
than station-keeping. Submersible simulators have not been employed to 
any great extent thus far, since the actual system is relatively inexpensive 
and the dynamics are usually very non-linear and difficult to model. The 
analysis of a submersible can be significantly more complex than the 
analysis of a conventional submarine or aircraft, since the presence of 
ancillary equipment such as manipulators, video devices, and tethers 
introduce extra cross-coupling terms usually absent in other, more 
symmetric, vehicle dynamic analyses. [Ref. 1] Additionally, all 
mathematical models include simplifying assumptions and errors in the 
model's hydrodynamic coefficients. 

Submersibles typically have a variety of complex dynamic interactions 
that can severely inhibit maneuverability and control performance. The 
goal of this thesis is to present an understanding of the coupling effects on 
straight lme motion stability in the horizontal plane using a hmnearized 
model, and the primary means of minimizing these effects. Development of 
the mathematical models for both the coupled and uncoupled maneuvering 
and roll equations of motion is presented in Chapter fl. Utilizing two 
different configurations of hydrodynamic coefficients, the degree of stability, 


regions of stability, and linear simulations for the coupled and uncoupled 


systems are presented in Chapter II]. A Hamming method nonlinear 
simulation [Ref. 2], which is similar to a fourth order Runge-Kutta 
integration technique, is conducted on both configurations to compare with 
the results obtained for the linearized models. Chapter IV develops a 
perturbation analysis to demonstrate the strong degree of dependence on 
the separation between the longitudinal centers of buoyancy and gravity to 
the solution of the linearized coupled equations of motion. Chapter V 
summarizes the results and provides recommendations for future 
submersible modelling research. Appendix A contains the computer 


programs utilized for the linear and nonlinear simulations. 


B. PARAMETER DEFINITION 
The values for the hydrodynamic derivatives and vehicle dimensions 
are from Smith, Crane, and Summey [Ref. 3], with the following exceptions: 
e Y, - the force in sway due to a change in yaw rate \ 
¢ N, - the moment due to a change in sway velocity. 
These two coefficients were modified to produce two different models that 
would have one eigenvalue change sign for a reasonable range of 


longitudinal and vertical centers of gravity. A comparison between the 


actual non-dimensional values and those used in this thesis is as follows: 





a= Ny 
Configuration ‘A’ = -3.500E-02 -1.484E-03 
Configuration ‘B’ _—-5.940E-02 -1.484E-03 
Actual SDV 2.970E-02 -7 420E-03 


The effects of changing these coefficients are illustrated in Figures 1 and 2 
on the following pages. Additionally, the analysis presented herein is 
conducted in dimensional form; hence the nominal forward longitudinal 


> 


velocity ‘U’ appears in the equations of motion. All calculations and 
simulations utilize a value of 5 ft/sec for ‘U’. The coordinate system 
convention is the standard body-fixed, right -hand orthogonal axis system 


employing the Euler angle approach. 


1. Variables 
meV, Z Distances along the body fixed principal axes. 
u, V, W Translational velocity components of model relative to 


fluid along body axes. 


Dp, a, 7 Rotational velocity components of model along body axes. 

mY, 2 Hydrodynamic force components along body axes. 

K, M,N Hydrodynamic moment components along body axes. 

y, 8, 0 Yaw, pitch, and roll angles; positive values following the 
right-hand rule. 

Xz, Yg, Ze, Center of gravity coordinates along principal axes. 

Ba, Yb, Zb Center of buoyancy coordinates along principal axes. 

Ixx, Lyy, zz Moments of inertia about the principal axes. 

Goce: stail Distances from body center along the longitudinal axis 
used in the crossflow force and moment integrals. 
Values are located within the nonlinear computer 
simulation program in Appendix A. 

h(x), b(x) Model width and height values used in the crossflow 


force and moment integrals. Values are located 
within the nonlinear computer simulation program 
in Appendix A. 
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II. STABILITY OF MOTION 


A. EQUATIONS OF MOTION 
The horizontal plane, nonlinear equations of motion for a submersible 
as developed by Smith, Crane, and Summey [Ref. 1] are shown below in 


Equations (2.1). 


Sway: m|[v + ur—wp+X¢(pq +f) — Ye(p? + r*) + Zg(qr- p)| = 
(2.1a) 
[Yop+ Yet + Ypqpq + Ygrar] +[Yyuv+ Yuavw+ Y5,u7r] + 


[Y,v + Ypup+ Y,ur+ Yyqvg + Ywpwp+ Ywrwr] + (W-B)cos@ sing - 


oe 9 9) (v+ xr) 
yy 
| [op »>h(x)(v + xr) + CDz b(x)(w -xq) ine dx 


Xtail c x) 


Yaw: Ti+ (ey ba : Beedle . I,,(pr+q) rt 
(2.1b) I,,(qr-p) + m[Xe(v+ur-wp) - Yg(u-vr+wq)} = 
[Npp + Nrr + Nogpq + Norar] + 
(Nyv + Npup + N,ur + Nyqgvq + Nwpwp +Nwrwy + 
(Nyuv + Nyww + N;,u6r] + (XgW-XbB)cos 6 sin d + (YgW- YbB)sin 6 - 
Tcpyh(x)(v 4x2 + Ope b(x)(w - x9) SEE (x) dx + u’Nprop 
X tail 


Nb Lpard, -1,,)+]x(pn-4) -ly.(q? - r7)- 
Ix2(pq +r) +m[Yg(w-uq+vp)-Zp(v + ur-wp)] = 
([Kpp+K;r+Kypgpq + Kogrqr] + 
[Kyv +K,up+K,ur+ Kygvq + Kypwp+Ky,wr] + 
[Kyuv+ Kywvw] + (Y,W - Y,B)cos 6 cos 6 - 


(ZgW - ZpB)cos 6 sin @ +u’Kprop 


B. SIMPLIFICATIONS 
In order to obtain the linearized equations of motion about a level flight 


path, the following simplifications were utilized: 


e The translational velocity (w) and acceleration (w ) in the z-direction 
are zero. 


e The rotational velocity (q) and acceleration (q) in the y-direction are 


zero. 


e The acceleration in the longitudinal direction (u ) is zero. 


e The cross-products of inertia are zero by virtue of a body-centered 
coordinate system. 


e The submersible is neutrally buoyant so B = W. 


e The longitudinal center of buoyancy (X,) and the vertical center of 
buoyancy (Z,) are located at the origin of the body-fixed coordinate 
system. 


¢ The lateral center of buoyancy (Yp) and the lateral center of gravity 
(Y,) are located at the origin of the body-fixed coordinate system. 


¢ Dynamic stability analysis is performed with all controls fixed; 
hence, all forces and moments due to control surfaces are zero. 


¢ The angle of pitch (8) is sufficiently small for sin(@) to equal zero. 


¢ From Smith, Crane, and Summey [Ref. 1] the propeller coefficients 
Kprop and Nprop are zero. 


C. COUPLED STABILITY EQUATIONS 
When the simplifying assumptions from the preceeding section are 


applied, the resulting linearized equations are: 


(2.2a) Sway: m{v +ur+X,(t)-Z,(p)] = Ys 
(2.2b) Yaw: L,t + mX,(v+ur) = Ng 
(2.2c) Roll: Ikxp - mZg(v+ur) = Kr 


where the force and moment representations are given by: 


(2.3a) Sway Force: Ye = Yyov + Yyv+ Yit+ Yrr+ Ypp+ Ypp 

(2.3b) Yaw Moment: Nr = Nyv +Nvyv+N;r+N,r+N,p+Npp+ (XgW- XpB)o 
(2.3c) Roll Moment: Ky = Kiv+Kyv+Kr+K,r+Kp+K,p+(Z,W-Z,B)o 

Equations 2.2 and 2.3 may be combined to form the state space 


representation: 


Ax = Bx (2.4a) 


where x=[p ov r]! (2.4b) 


Ixx-Kp 0 -(Kv+ MZ) -Kr 


-(Yp+MZg) 0 M- Yv MX¢- Yr 


-Np 0 MX-Nv Izz - Nr 


KpU ZbB-ZgW KvU U(MZg+ Kr) 


1 0 0 ) 
B = 2.4d 
YpU 0 YvU U(Yr-M) ee 


NpU XgW-XbB NvU U(Nr-MXzg) 

The stability of the coupled, linearized system depends on the location of the 
four eigenvalues of det(B - AA) = 0, which has a characteristic equation of 
the form 

Al4+BA34+CA2+DA+E=0, (2.4e) 
where the coefficients are complicated permutations of the elements in 
matrices A and B. The values for A, B, C, D, and E are given in Equations 
(2.4f - 2.4}) in terms of lower case letters that represent entries in matrices A 


and B; they are explicitly defined in Table 1. 


A = a-(j-u - |-r) + d-(u-h + o-l) - f-(r-h + 0-j) (2.4f) 
B= e-(u-h + 0-1) - d-(u-i - w-h + o-m - p-l) - a-G-w + k-u - I-x - r-m) 

- b-(j-u - l-r) + f-(r-i - x-h + o-k - p-j) - g-(r-h + 0-)) (2.4¢) 
C = a-(k-w - m-x) + b-(Gj-w + k-u - ]-x - r-m) + g-(r-i - x-h + o-k - pj) 

+ ¢c-(j-u - l-r) - d-(q-] - w-i + m-p) - e-(u-i - w-h + o-m - p-l) 


+ f-(q-:j - x-i + p-k) (2.4h) 


D= g-(q:j - x-i + p-k) - f-q-k + d-q-m - e-(q-] - w-i + m-p) - b-(k-w - m-x) 
- ¢c-(j-w + k-u - 1-x - rem) (2.41) 


E= c-(k-w - m-x) + e-q-m - g-q-k (2.4}) 


TABLE 1. COEFFICIENT DESCRIPTOR VALUES 





c = ZgW- ZB MZg + Ky 
g = MZgU MZg + Yp 
k= MXg - Yr 


m = U(Yr-M) NpU XbB - XgW 


r = MxXg- Nv U(Nr-MXg) x NvU 


D. UNCOUPLED STABILITY EQUATIONS 
Using matrices A and B from the preceeding section (Equations 2.4c and 


2.4d), uncoupling the steering and roll equations is straightforward: 


STEERING EQUATIONS 
M- Yv 2 E : eS U(Yr-M) I] 


MXg-Nv Ix-Nr |[r NvU U(Nr-MXg)] Lr (2.5a) 
ROLL EQUATIONS 
ors : ia 7 ire 4 P| 
0 1 eee 0 Ld (2.5b) 


10 


The characteristic equations for the uncoupled steering and roll equations 
are much simpler than that for the coupled system, and are given as 
Equations (2.5c and 2.5d) in terms of the hydrodynamic parameters rather 


than the descriptors used in the previous section. 


The steering characteristic equation form is: Ap + BrA + Cy = 0, where 


A, = (M-Yy lw -Nz) - (MX, - Y; (MX, -Ny) 


BL = -[dz-N; (YU)+(M- Yy (Ny - MX, (U)+ (2.5c) 
(Yr -M)(Ny -MXg (U) + CY; -MXg (NvU)] 


Cee. (Neen) - (NU (Y= 


The roll characteristic equation form is: eae + Bra + Cr = 0, where 


Ar = kx-K, 
Br = -K,U (2.5d) 
Cr = Z,W 
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III. RESULTS AND DISCUSSION 


A. DEGREE OF STABILITY 


The effects of changing the longitudinal and vertical centers of gravity 
(while keeping the center of buoyancy fixed at the vehicle center) on the 
degree of stability for configuration ‘A’ is illustrated in Figure 3. Degree of 
stability as utilized in this thesis is defined as the maximum real value of all 
characteristic roots, with negative values indicating a stable situation. The 
degree of stability for the uncoupled system is represented by the dashed line. 
It can be seen that the critical value of X, for which motions become unstable 
is clearly a function of the metacentric height Z,, whereas the uncoupled 
system predicts a constant value of X,. Figure 4 displays the variation of the 
imaginary part of the root value for configuration ‘A’. 

Figures 5 and 6 are analogous to Figures 3 and 4 for configuration ‘B’; 
they show the degree of stability for varying X, and Z, values. For this 
configuration, the degree of stability has a stronger dependence on the 
location of X, and the value of Z,. For almost all positive values of X,, the 
complex conjugate roots are increasing in value and eventually becoming 
positive; this indicates an oscillatory response that diverges when the degree 


of stability is positive. 


B. REGIONS OF STABILITY 

Figure 7 shows the region of stability for configuration ‘A’, with the 
uncoupled system represented by a dashed line. The uncoupled system 
predicts stability for all values of X, greater than 0.18 ft, while the coupled 


system predicts an additional region enclosed by the triangular area to the 
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left of the dashed line. Figure 8 is analogous to Figure 7 for configuration ‘B’. 
The large discrepancies between the predicted regions of stability in this case 
occur for Z, values less than 0.045 ft. For small values of Z,, there are 
corresponding small regions of X, where stability is predicted by the coupled 
equations but not by the uncoupled equations. The root values in this region 
are complex conjugates with very small real parts. Figures 9(a) and 9(b) 
illustrate the effect of co-locating Xp and X,. This results in a degree of 
stability for the coupled system that is nearly identical to the uncoupled 
system. 
C. INTERPRETATION OF RESULTS 

It may be shown by applying Routh’s criterion [Ref. 4: pp. 211-218] that 
for a fourth order equation of the form AA* + BA® + CA? + DA + E to have 
roots with all negative real parts the following must apply: 

i.) BCD ABs aBeet 
ii.) E > 0. 

If the quantity ‘E’ is less than zero, the system will become unstable and the 
resulting motion will be a simple divergence. If, however, the value of the 
quantity BCD - AD? - EB” is negative, the resulting instability will result in 
an oscillatory motion due to the presence of complex conjugate roots with 
positive real parts. 

For the coupled system of equations, the condition E = 0 yields: 

Ze = X_-[KAY-M)/(YN, - NAYM))]. 

while the uncoupled system of equations reduces to a constant term 
expression for Xg: 


X, = LOMy,)/(YN, - NUY-M))]. 
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This explains the differences in the regions of stability illustrated in Figure 7, 
since the above equations show a linear relationship between Z, and X, for 
the coupled equations and a constant value for X, for the uncoupled 
equations. When the value for X, coincides with the value for Xp, the 
constant term ‘E’ in the coupled equations is reduced to that of the constant 
associated with the uncoupled equations, and the resulting predicted degree 
of stability no longer depends on the value of Z,. Substituting the coefficient 


’ 


values for ‘E’ serves to clearly illustrate the reduction: 
E = (Z,W)L(YU2)(Ny-MX,) - (NWU2¥,-M)] + 
(K\U*)(Y,-M\(XpB-XgW) - (MZ,U)YU)(XpB-X,W). 


When Xj = X, and B = W for the neutrally buoyant case, the second and third 
terms are reduced to zero. When ‘E’ is then set equal to zero (the condition 
for determining where the real roots change from negative to positive), the 
dependence on the value Z, is removed and the expression for X, is the 
constant given for the uncoupled equations. This reduction is the explanation 
for the appearance of Figures 9(a) and 9(b). When the longitudinal centers of 
buoyancy and gravity coincide for a neutrally buoyant vehicle, the degree of 
stability for the coupled system of equations covering all metacentric height 
values is equivalent to the uncoupled system of equations. 

A simple reduction of the equation resulting from Routh’s criterion to 
determine when a pair of complex conjugate roots crosses the zero axis is not 
easily accomplished. Figure 10 is presented as confirmation that the 
locations of X, for which BCD - AD? - EB” = 0 matches the locations given 


graphically in Figure 5. 
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D. LINEAR SIMULATIONS 

Figures 11 through 14 present comparisons of the coupled and uncoupled 
system responses for configuration ‘A’. For the stable cases X, = 0.40 ft and 
for the unstable cases X, = -0.20 ft, while Z, is 0.20 ft for both. The unstable 
coupled case (Figure 11) illustrates a simple divergence for both angle of drift 
and roll angle. This is expected since the roots are not complex conjugates. 
The unstable, uncoupled case accurately predicts the divergence in angle of 
drift, but the roll response is predicted to be stable. This may be explained by 
examining the uncoupled equation of motion in roll: 

O° - o' (KpU)/(lx - Kp) + 6 (ZgW)/(lx - Ky) = 0. 
Substituting values for configuration ‘A’ yields: 
o” + o° (1.475) + $ (0.720) = 0. 
The solution of this equation results in a natural frequency of 0.85 rad/sec 
and a damping factor of 0.869. This is an underdamped case, since both roots 
have negative real parts and are complex conjugates. Substituting values for 
configuration ‘B’ yields: 
o” + 0’ (1.475) + (0.144) =0. 

The solution for this case results in a natural frequency of 0.380 rad/sec and a 
damping factor of 1.941, which represents an overdamped situation. This 
also demonstrates that a vehicle’s natural frequency in roll may be increased 
by increasing the metacentric height. The effects of the underdamping may 
be seen in Figures 12 through 14. 

Figures 15 through 20 provide a comparison between the coupled and 
uncoupled equations for configuration ‘B’. The effects of the overdamping are 


evident in figures 15 through 17. This set of figures demonstrate the 
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magnitude of the discrepancies when the coupling effects are not considered. 
The comparison below summarizes the results of the figures representing 
simulations for configuration ‘B’, where ‘C’ stands for coupled and ‘UC’ for 


uncoupled. 


C UC C UC C UC 


X, 0.20 0.20 1.00 1.00 T50 1.50 
Roll Stable Y ny N Y N NE 
Drift Stable ey. N N N N bg 


As was seen before, the effects of the coupling on the system results in 
complex conjugate eigenvalues with positive real parts. Therefore, the larger 
values of X, result in increasingly divergent oscillations instead of the 
stability predicted by the uncoupled system. 

Figure 21 is a three-dimensional presentation of the roll amplitude vs 
time for configuration ‘B’ as X,g varies from 0.15 ft to 1.50 ft. This mesh 
capability of MATLAB allows a comparative view of several solutions, and 


the behavior of the roll response for the coupled system is easier to discern. 


E. NON-LINEAR SIMULATIONS 

In order to provide a measure of the accuracy of the results obtained 
utilizing the linearized equations of motion, a simulation program for the 
non-linear equations of motion was developed using Hamming’s method 
[Ref. 2]. Hamming’s method utilizes a Milne predictor and incorporates a 
modifier step prior to the correction step. The primary advantage of using 


Hamming’s method is that only two derivative function evaluations are 
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required per step rather than the four or more normally required by other 
popular methods. The local error is of the same order of magnitude (h‘) as a 
more time consuming process such as a fourth order Runge-Kutta, but the 
reduction in function evaluations results in a faster simulation. The formula 
is presented below, and may also be found in the nonlinear computer 


program simulation in Appendix A. 


HAMMING’S METHOD 
yli+1)predicted = y(i-3) + (4h/3)[2f0i) - £G-1) + 2fG-2)] 
Ylit+] modified = Ylit1)predicted - (112/121) YG) predicted - YCi)correctea] 
yi+1)corrected = (1/8)[9yCi) - yG-2) + Sh{fi+ 1) modified + 2Mi) - fi-1)}] 
yG+1) = ylitD corrected + (9/121 [yli+ predicted - ylit]) corrected] 
The first four values must be determined by another method; the Euler 
linear solution with a small step size ‘h’ proved sufficient. 

Figure 22 represents the simulation for the unstable representation of 
configuration ‘A’. Rather than the exponentially increasing roll angle and the 
-90 angle of drift computed with the linear simulation, the nonlinear solution 
predicts an angle of drift that reaches -15 degrees, and then slowly diverges. 
The roll angle reaches a steady state value of approximately three degrees. 
Figures 23 and 24 are the nonlinear simulation results for stable 
configurations of ‘A’ and ‘B’, respectively. They are nearly identical to the 
results obtained using the linear simulation and displayed as Figures 12 and 
18. 

Figure 25 is the simulation for an unstable configuration ‘B’. Quite 
notable are the steady state roll angle and angle of drift after approximately 
250 seconds rather than the exponentially increasing divergence apparent in 
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the corresponding linear simulation. Figure 26 provides a comparison 
between the linear and nonlinear solution for angle of roll. The enlarged 
view of the steady state region provided in Figure 27 better illustrates the 
constant limits of oscillation. By plotting the roll angle versus the angle of 
drift, an elliptical trajectory is apparent in Figure 28. This result is similar 
to the results obtained by Schmidt and Wright in their analysis of high 
performance aircraft wing-rock [Ref. 5]. They postulate that a possible 
explanation for the limit cycle is the inertial coupling between a stable 
longitudinal and an unstable lateral mode. Similar results in the work are 
attributed to dynamic as well as hydrodynamic coupling. Figures 25 and 28 
indicate that the nonlinear interactions for the unstable conditions of 
configuration ‘B’ provide a significant amount of damping to the rolling 
motion. The limiting of the rolling motion accomplished by including the 
nonlinear terms then serves to limit the buildup of the angle of drift and the 
result is an eventual stable limit cycle. 

The ocean environment can be expected to introduce many combinations 
of disturbance forces, which may or may not be periodic. A preferred method 
for simulating many disturbances such as sensor noise, ocean current 
fluctuation, and vehicle acceleration fluctuations is to model them as ‘white 
noise. Figure 29 is presented to demonstrate the effects of including 
constant, zero-mean disturbances in sway and yaw on the marginally stable 
system of configuration ‘B. The disturbances are developed using MATLAB’s 
random number genéri:tor with a uniform dij:tribution. Che random numbers 
are then scaled to simulate values that may be < «pected. The variance of the 
sway and yaw acceleration disturbances, respectively, for Figure 29 are 0.003 


(ft/s”)” and 0.002 (rad/s”). The resulting simulation bears little resemblance 
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to Figure 24, although the initial conditions and values for X, and Z, are the 
same. With even these relatively minor disturbances acting on the system of 
configuration ‘B’, a rather large, non-symmetric oscillation in both roll and 
angle of drift is evident. Although the system is still stable, with the mean of 
both the roll angle and angle of drift equalling zero, a limit cycle similar to 
that of the unstable configuration (Figure 25) has developed. As the angle of 
drift fluctuates between positive three degrees and negative four degrees, the 
angle of roll varies between positive and negative two degrees. Increasing 
the scaling (which increases the variance) for the disturbances would be 
expected to increase the fluctuations until stability is lost. Similarly, it is 
true that the small disturbances acting in Figure 29 have a greatly reduced 


effect on the system of configuration ‘A’, which has greater stability. 
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IV. PERTURBATION ANALYSIS 


A. BACKGROUND 

The previous section demonstrates that knowledge of the separation 
between longitudinal centers of buoyancy and gravity is critical in 
determining system stability. If this quantity is the dominant factor in a 
stability analysis, then an approximation for the degree of stability may be 
developed by applying perturbation theory. Although perturbation results 
are only approximations, their advantage over numerical methods is in 
illustrating the degree to which a solution depends on the variable(s) 


involved. From the fundamental theorem of perturbation theory [Ref. 6], we 
seek a solution to the characteristic equation Ak* + BA’ + Ci? + DA + ESO 


the form: 


Me) = y an e” (4,1) 
n=0 


where € is the variable of interest and the solution approximates the 
numerical solution in the region where € is smaill. The constant coefficients 
(ao, a1, ..-, An) are all independent of €, and are all equal to zero for small e€. 
B. PERTURBATION FORMULATION 

Recalling Equations 2.4f - 2.4j, which defimed the coefficients of the 
quartic characteristic equation, it would obviously be desirable to simplify the 


equations as much as possible prior to formulating the perturbation 


approximation. The variable of interest will be defined as: 


e=X,- Xp (4.2) 
where the nominal value for X, is zero and the perturbation will] be 


performed around ec = 0. By comparing Figures 30 and 31 it is clear that by 


allowing the hydrodynamic coefficients Kv, Ky, Y;, Np, Yp, and N to equal 


zero a very good approximate solution to Equation (2.4e) is obtained for both 
configurations ‘A’ and ‘B’. This simplification reduces the descriptor 
coefficients ‘f, ‘i’, ‘o’, and ‘p’ from Table 1 to zero, which simplifies the 


coefficients of the coupled equations of motion to Equations (4.3a) - (4.3e) 


below. 


A = (xx - Kp)[(M- Yy)(zz -N+)-(MXg- Y;)(MXg-Ny)} + 
(MZ, )*(Iz2 - N;) (4.3a) 


B = (KyU)(MZg)(Iz2 - Ny) + (MZ,)?(N; -MX,)(U)- 
(MZ,)(MZg +Kr)(MXg - Ny )(U)- Tux - Kp X(M- Yy (Ny - MX g)(U) + 
(Y,U)(Iez - Nz) - (NVU)(MX g - ¥;)- (MX - Ny (Yr - M)(U)]- 
(Kp U){(M - Yy (Ize - Nr )- (MXg - Y; (MX g- Ny) (4.3b) 


C = (Ixx - Kp)[(YyU?)(Nr -MXg)-(NvU*)(¥,- My + 
(ZgW(M - Yu )(Izz - N;)-(MXg - ¥-)(MX, -Ny)] + 
(K,U)E(M - Ys )(Nr- MX, )(U) + (¥, U)(Ie2 -N;) - 
(N\U)\(MXg - Y+) - (MX_-Ny)(¥; -M)(U)] + 
(MZ, )[(X,W)(MXg -Y;)-(N,U?)(MZ, + K,) +(KyU?)(N; -MXQ)] 


(4.3c) 


D = -(Z,W)I(M- Yy)(N, -MX,)(U) +(YyU)Uz2 -N¢)- 
(N vVU)(MX ¢ - Y;)-(MXg - Nv)(¥r-M)(U)] - 
(KpU®)[(Y, (N; - MXg)-(Ny)(¥;-M)] + (KyU)(XgW)(MXg - Y;)- 
(MZ, )(X,W)(U)(Y; - M) - (MZ, +K-)(U)(X, W)(M - Yy ) (4.3d) 


G3] 
I 


(ZgW)(Y,U?)(N, -MX,)-(NyU?)(¥,-M)] - (KyX,WU?)(Y, -M) + 
(YyX gWU?)(MZe + Ky) (4.3e) 


Recalling the definitions of the coefficients for the uncoupled systems: 
AL = (M - Yy)(z2z -N;) . (MX¢g Y; (MX g -Nv) 


Br, = - ee -N;)(¥vU) +(M- Yy (Nr - MXg)(U) + 
OG = M\(Ny -MX,)(U) + (Y; B MX, \N,U)] 


Cy = (YyU?)(Nr-MXg) - (NvU?)\(Y; -M) 


AR = Ixx-K, Br = -K,U Cr =Z,W 


Substituting in Equations (4.38a) - (4.3e) yields: 


A = ApAy + (MZg)*(Izz -N;) (4.4a) 
B = AgBy + BrAy + (MZ,z)*(N, -MX,)(U) + (4.4) 

(MZez Y(Iz2 - Ni)(KyU)- (MZ, + Ky (MX, - Ny )(U)] 
C = ArpC, + CrAz + BRB + (MZ, )I(XgW)(MX,g - Y;) + 

(K U2)(N; - MX,) -(NyU2)(MZ, + K,)] ae 

D = CrBy + CLBR + (KyU)(X,W)(MX, - Y;) - 

(MZ, + Ky)(XgW)U)(M- Yy) - (MZg)(XgW\U)(Y; -M) (4.44) 
E = CrCy + (YvU*)(X,W)(MZe + Kr) - (KvXgWU?)(Yr -M) ae 
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In order to further simplify these coefficients, terms of order (Zp) and (Xa) 


will be neglected; the effects of doing this are small and may be seen in 


Figures 32 and 33. This reduces the coefficients to: 


A = ApAz (4.5a) 
B = ApBL + BRA + a (X,) + K1 (4.5b) 
C = ApCy + CpAy + BRB + B (X,) + K2 (4.5c) 
D = CRB, + CLBp + (X,) (4.5d) 
E = CrCy + 8 (XQ) (4.5e) 
where: o = -M°Z,K,U (4.6a) 
B = -(MZg)[(WY:) + (MKvU*)] (4.6b) 
y = (WUYMZg (We -Yr)+ Kr¥e-KrM-KvYr] (4.6c) 
& = (WU”)[(Y,)(MZ,+K;) - (Ky)(Y¥,-M)] (4.6e) 
K1 = (MZg){(Izz- Nr )(KvU) + (NvyKrU)] (4.6f) 
K2 = (MZ,U*)(N,K, - NyK,) (4.6g) 


Carrying through the computations results in the following expressions: 
(ARAL) ao’ + (ARBL+BRALtK1) ao? + (ApCL+CrAL+BRBy+K2) ao” + 
(CpB,+CLBpR) an + CRC, = O (4.7a) 


~(aao° a Bao? + Yao + 8) 





(4.7b) 


“1 © (GA)ao? + 3(B- aX, aoe + 2(C-BX,)a0 + (D- yg) 
Es -a17(C- BXz) - a;7ao(6A +B-aX,) - K laoa?’ - 30a02a1 - 2Baoa1- a1 
2 oe nen ee EEE 


(4A )ao® + 3(B -0Xg)ao” + 2C-BXg)ao +(D- yXg) 
(4.7c) 
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The characteristic equation AA* + Bh 4 Ch* 4 eee) may now be 


expressed in terms of a second order perturbation which depends on a single 


variable € as: 


Me) = ag + a, € + a e7 (4.8) 


Figures 34 through 37 depict the results for both first and second order 


perturbation analysis about X, close to zero. Results for both configurations 
‘A’ and ‘B’ demonstrate that this method for representing the degree of 
stability in the vicinity of a nominal (Xg - Xp) is quite satisfactory, provided 


the roots have no complex conjugates. In the region where the roots become 
complex, the problem is no longer a regular perturbation, and methods to 


solve the singular perturbation must be developed. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


When considering an inherently stable situation for an analysis of sway, 
yaw, and roll response, the linear simulations compare favorably with 
the nonlinear simulations. For a marginally stable or unstable system, 
only the nonlinear simulation can predict the existence and extent of any 
non-zero steady states or limit cycles that may occur. 

Obviously, in the design of a submersible, careful attention must be paid 
to selection of both the metacentric height and the separation between the 
longitudinal centers of gravity and buoyancy. Instability associated with 
the dynamic coupling effects can be minimized by increasing the 
metacentric height. The uncoupled system stability predictions are not 
reliable when there is separation between X, and Xp. 

It can be concluded from the analysis that the coupling of roll into sway 
and yaw for the linearized equations of motion is not very significant 
when Xz equals Xp. 

Recommendations for future modelling research include an expansion 
to incorporate coupling effects between the vertical and horizontal planes. 
It is well known that the coupling between translation and attitude is a 
severe limitation in the design of submersibles. 

The effects of varying hydrodynamic derivatives and initial conditions, as 
well as incorporating disturbances in the nonlinear simulations deserve 
more attention. Additionally, the complexity of the nonlinear system 
resulting from including the vertical plane is greatly increased. More 
studies of the nonlinear system resulting from this union are also 


required. 


APPENDIX A 
N-LINE IMULATION P 


AK Ak 2 2k 2 2 2 2 2K fe 2k 2 2k OK 2 2k 9K 6 2K 2 kK 2K ok OK 2K 9K OK 2 OK OK 6 KOK 2 OK OK 2 KOK OK OK OK 2 2 OK OK 2 OK OK 2 KK OK KOK KK KKK KK KK 


% THIS PROGRAM USES HAMMING’S METHOD TO SOLVE A SYSTEM 
% OF NONLINEAR EQUATIONS. SIMILAR TO A FOURTH ORDER 

% RUNGE-KUTTA IN TERMS OF LOCAL ERROR, IT REQUIRES ONLY 
% TWO FUNCTION CALCULATIONS PER STEP VICE FOUR. 


function[v,phi] = nonlinsim_hamming(U0,Xg,Zg) 
% VEHICLE PARAMETERS FOLLOW 
W = 12000.0; Ixx = 1760.0; Iyy = 9450.0; Izz=10700.0; L = 17.425; 
RHO=1.94; G=32.2; M = W/G; BUOY = W; Zb=0.0; Xb = 0.0; 
7Zi=0.5* RAO 2: 
% SWAY HYDRODYNAMIC COEFFICIENTS 


Ypdot = 1.270E-04*ZZ*L“2; Yrdot = 1.240E-03*ZZ*L‘2; 
Ypq = 4.125E-03*ZZ*L‘2; Yqr = -6.510E-03*ZZ*L“2; 
Yvdot =-5.550E-02*ZZ*L; Yp = 3.055E-03*ZZ*L; 
Yvq = 2.360E-02*ZZ*L; Ywp = = 2.350E-01*ZZ*L; 
Ywr =-1.880E-02*ZZ*L; Yv = -9.310E-02*ZZ; 
Yvw = 6.840E-02*2ZZ,; 


% NOMINAL VALUE FOR Yr = 2.970E-02*ZZ*L 
% CONFIGURATION ‘A’: Yr = -3.500E-02*ZZ*L 
% CONFIGURATION ‘B’: Yr = -5.940E-02*ZZ*L 


CC = input(‘enter model configuration choice; either 1 for A or 2 for B’) 


if CC <2 

Yr = -3.500e-02*ZZ*L; 
else 

Yr = -5.940e-02*ZZ*L:; 
end 


% ROLL HYDRODYNAMIC COEFFICIENTS 


Kpdot = -1.01E-03*ZZ*L‘3; Krdot = -3.37H-05*ZZ*L‘3; 
Kpq = =-6.93E-05*ZZ*L‘3; Kqr = 1.68E-02*ZZ*L‘3; 
Kvdot = 1.27E-04*ZZ*L‘2; Kp = -1.10E-02*ZZ*L“2; 
Kr = -8.41E-04*ZZ*L“2; Kvq =-5.115E-03*ZZ*L“2; 
Kwp = =-1.27E-04*ZZ*L‘2; Kwr = 1.89E-02*ZZ*L“2; 
Kv = 3.055E-03*ZZ*L; Kvw =-1.87E-01*ZZ*L; 


&2 


% YAW HYDRODYNAMIC COEFFICIENTS 


Npdot = -3.370E-05*ZZ*L‘3; Nrdot =-3.400E-03*ZZ*LA3: 
Npq = -2.110E-02*ZZ*L‘3; Nar = 2.750E-03*ZZ*L3; 
Nvdot = 1.240E-03*ZZ*L/2; Np = -8.405E-04*ZZ*L/2; 
Nr = -1.640E-02*ZZ*L/2; Nvq  =-9.990E-03*ZZ*L/2; 
Nwp =-1.750E-02*ZZ*L/2: Nwr = 7.350E-03*ZZ*L2; 
Nv = -1.484E-02*ZZ*L; Nvw = -2.670E-02*ZZ*L; 


% THE FOLLOWING ARE USED FOR EVALUATING THE INTEGRALS 
% IN THE SWAY AND YAW EQUATIONS OF MOTION. ‘X’ IS THE 

% DISTANCE IN FEET ALONG THE LONGITUDINAL AXIS AND ‘HH’ 
% 1S THE VEHICLE HEIGHT. ALL VALUES ARE IN FEET. 


X(1)=-105.9/12; X(2)=-104.3/12; X(3)=-99.3/12; 
X(4)=-94.3/12; X(5)=-87.3/12; X(6)=-76.8/12; 
X(7)=-66.3/12; X(8)=-55.8/12; X@)=72.7/12; 
X(10)=79.2/12; X(11)=83.2/12; X(12)=87.2/12; 
X(13)=91.2/12; X(14)=95.2/12; X(15)=99.2/12; 
X(16)=101.2/12; X(17)=102. 1/12; X(18)=103.2/12; 
HH(1)=0.0/12; HH(2)=2.28/12; HH(3)=8.24/12; 
HH(4)=13.96/12; HH(5)=19.76/12; HH(6)=25.1/12; 
HH(7)=29.36/12; HH(8)=31.85/12; HH(9)=31.85/12; 
HH(10)=30.00/12; HH(11)=27.84/12; HH(12)=25.12/12; 
HH(13)=21.44/12; HH(14)=17.12/12; HH(15)=12.0/12; 
HH(16)=9.12/12; HH(17)=6.72/12: HH(18)=0.00/12; 


% THE INITIAL CONDITIONS ARE SET PRIOR TO INTEGRATION 
% BY CALLING A LINEAR INTEGRATION PROGRAM NAMED 

% ‘EULER’. THE ERRORS ASSOCIATED WITH USING LINEAR 

% SOLUTIONS FOR THE FIRST FIVE TIME INTERVAL STEPS ARE 
% SMALL AND DO NOT AFFECT THE NONLINEAR SOLUTIONS. 


[p,pdot,v,vdot,r,rdot,phi,phidot] = euler(dt); 


p(1:5)= p; 
pdot(1:5)= pdot; 
pdotmod(5)=pdot(5); 


v(1:5)= v; 
vdot(1:5)= vdot; 
vdotmod(5)=vdot(5); 


Heo =. 
rdot(1:5)= rdot: 
rdotmod(5)=rdot(5); 


phi(1:5)= phi; 
phidot(1:5)= pludot; 
phidotmod(5)=phidot(5); 


% THE TIME INTERVAL IS ‘dt’, THE FINAL TIME IS ’tfinal’, AND A 
% VALUE CALLED ‘stopnumber’ IS SET TO ALLOW ACCESS TO THE 
% DATA AFTER THIS NUMBER OF ITERATIONS. THE VALUE FOR 
%o ‘stopnumber MAY BE CHANGED WHILE THE KEYBOARD IS 

% ACTIVE TO ALLOW SUBSEQUENT PROGRAM INTERACTION. 

% RETURN TO THE PROGRAM IS ACHIEVED BY TYPING ‘return’ 
% FOLLOWED BY THE ‘enter’ KEY. 


dt = input(‘enter the time interval step size’) 
tfinal = input(‘enter the final time’) 
stopnumber = input(‘enter the value for stopnumber’) 


% THIS SECTION ALLOWS FOR RANDOM DISTURBANCES IN SWAY 
% AND YAW 
rand(‘uniform’) 
vdotdist = rand(1,(tfinal/dt)+1); 
vdotdist = 0.05*(vdotdist - mean(vdotdist)); 
rdotdist = rand(1,(tfinal/dt)+1); 
rdotdist = 0.04*(rdotdist - mean(rdotdist)); 


% THIS SECTION PROVIDES FOR CONTINUATION OF SOLUTIONS. 
% IF THE MATRICES BECOME VERY LARGE, IT IS RECOMMENDED 
% THAT THE CURRENT VALUES BE SAVED AND A NEW 

% SIMULATION BEGUN AS MATLAB SLOWS NOTICEABLY WITH 
% LARGE MATRICES. 


% “Vcorri(5)= rcorr(5)= } 
% pcorr(5)= phicorr(5)= ; 
% vpred(5)= rpred(5)= ; 
% ppred(5)= phipred(5)= ; 
% pdotmod(5)=  ; rdotmod(5)= ; 


% COMMENCE INTEGRATION 


for ) = 6:(tfinal/dt) + 1; 
% THIS SECTION USES A HAMMING PREDICTOR-CORRECTOR TO 


OBTAIN VALUES 


vpred(j)=v(j-4)+(4*dt/3)*(2*vdot(j-1)-vdot(j-2)+2*vdot(j-3)); 
rpred(j)=r(j-4)+(4*dt/3)*(2*rdot(j-1)-rdot(j-2)+2*rdot(j-3)); 
ppred(j)=p(j-4)+(4*dt/3)*(2*pdot(j-1)-pdot(j-2)+2*pdot(j-3)); 
phipred(j)=phi(j-4)+(4*dt/3)*(2*pG-1)-pG-2)+2* p(j-3)); 


Hy< 7 
vmod(j) = 
rmod(j) = 
pmod(j)= 
phimod(j)= 

else 
vmod(j)= 
rmod(j)= 
pmod(j)= 
phimod(j)= 
end 


vpred(j); 
rpred(j); 
ppred(j); 
phipred(j); 


vpred(j)-(112/121)*(vpred(j-1)-vcorrG-1)); 
rpred(j)-(112/121)*(rpred(j-1)-rcorr(Gj-1)); 
ppred(j)-(112/121)*(ppred(j-1)-pcorr(j-1)); 
phipred(@j)-(112/121)*(phipred(j-1)-phicorr(j-1)); 


tmp1mod=vmod(j); tmp2mod=rmod(j); 
[SWAYMOD, YAWMOD]=crossflowintmod(tmp1mod,tmp2mod,X,HH); 


vdotmod(j)= 


vcorr(j)= 


rdotmod(j)= 


rcorr(j)= 


pdotmod(j)= 


(1/((M-Yvdot))*((Ypdot+M*Zg)*pdotmod(j-1)+... 
Yp*U0*pmod(j)+(Yrdot-M*Xg)*rdotmod(j-1)+... 
Yv*U0*vmod(j)+(Yr-M)*U0*rmod(j)+SWAYMOD); 


0.125*(9* v(j-1)-vG-3)+3*dt*(vdotmod(j)+2*vdot(j-1)-... 
vdot(j-2))); 


(1/(Izz-Nrdot))*(Npdot*pdotmod(-1)+... 
(Nvdot-M*Xg)*vdotmod(j)+Np* U0*pmod(j)+... 
(Nr-M*Xg)*U0*rmod(¥j)+Nv*U0*vmod(j)+... 
Xg*W*sin(phimod(j))+YAWMOD); 


0.125*(9*rG-1)-rG-3)+3*dt*(rdotmod({j)+2*rdotG-1)-... 
rdot(j-2))); 


(1/(Iixx-Kpdot))*((Kvdot+M*Zg)*vdotmod(j)+... 
Krdot*rdotmod(j)+Kp*U0*pmod(})+... 
(Kr+M*Zg)*U0*rmod(j)+Kv*U0*vmod()-... 
Zg*W*sin(phimod(j))); 


peorr(j)= 0.125*(9*p(j-1)-pG-3)+3*dt*(pdotmod(j)+2*pdot(j-1)-... 
pdot(j-2))); 

phicorr(j)= 0.125*(9*phi(j-1)-phig-3)+3*dt*(pmod(Gj)+2*pGj-1... 
p(j-2))); 

v= veorr(j)+(9/121)*(vpred(j)-vcorr(j)); 


rQ)= 
pQ= 


rcorr(Gj)+(9/121)*(rpred(j)-rcorr(j)); 
pceorr(j)+(9/121)*(ppred(j)-pcorr(j)); 
phicorr(j)+(9/121)*(phipredG)-phicorr(j)); 


betaQ) = -atan(v(j)/5)*180/pi; 


% THIS SECTION HAS THE NON-LINEAR EQUATIONS OF MOTION. 
% REMEMBER TO REMOVE THE DISTURBANCES FROM THE SWAY 
% AND YAW DERIVATIVES IF THEY ARE NOT DESIRED. 


tmpl=v(j); tmp2=r(j); 
[SWAY,YAW] = crossflowint(tmp1,tmp2,X,HH); 


vdotG)= ((1M-Yvdot))*((Ypdot+M*Zg)*pdotG-1)+Yp*U0*p(j) +... 
(Yrdot-M*Xg)*rdotG-1)+Yv*U0*vG)+(Yr-M )* U0* r(j)+... 
SWAY))+vdotdist(j); 

rdot(j= ((1/Izz-Nrdot))*(Npdot*pdotG-1)+(Nvdot-M*Xg)*vdotG}+... 
Np*U0*pGj)+(Nr-M*Xg)* U0*r(j)+Nv* U0*v(j)+... 
Xg*W*sin(phi(j))+YAW))+rdotdist(j); 

pdot(j= (1/(Ixx-Kpdot))*((Kvdot+M*Zg)*vdot(j)+Krdot*rdot(j)+... 
Kp*U0*pG)+(Kr+M*Zg)*U0*rGj)+Kv*U0*v(j)-... 
Zg*W*sin(phi(j))); 

if } == stopnumber 

keyboard 

end 

end 

return 
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% THIS PROGRAM IS THE LINEAR EQUATION OF MOTION 
% INTEGRATION EMPLOYING A SIMPLE EULER METHOD TO 
% OBTAIN STARTING VALUES FOR THE NONLINEAR SIMULATION. 


function[p,pdot,v,vdot,r,rdot,phi,phidot] = euler(dt) 


% THE MATRIX ‘A’ CORRESPONDS TO THE DAMPING MATRIX AND 
% MATRIX ‘B’ CORRESPONDS TO THE MASS MATRIX. 


A(1,1) = Kp*U0; 
A(1,3) = Kv*U0; 


A(1,2) = -((Zg*W)-(Zb*B)); 
A(1,4) = (Kr+(M*Zg))*U0; 


P(Z,1) = 1.0; A(2,2) = 0.0; 
A(2,3) = 0.0; A(2,4) = 0.0; 
A(3,1) = Yp*U0; A(3,2) = 0.0; 


A(3,3) = Yv*U0; 
A(4,1) = Np*U0; 
A(4,3) = Nv*U0; 


B(1,1) = Ixx-Kpdot; 

B(1,3) = -(Kvdot+(M*Zg)); 
B(2,1) = 0.0; 

B(2,3) = 0.0; 

B(3,1) = -(Ypdot+(M*Zg)); 
B(3,3) = M-Yvdot; 

B(4,1) = -Npdot; 

B(4,3) = (M*Xg)-Nvdot; 


C = [inv(B)*A); 
% INITIAL CONDITIONS 


A(3,4) = (Yr-M)*U0; 
A(4,2) = ((Xg*W)-(Xb*B)); 
A(4,4) = (Nr-(M*Xg))*U0; 


Bo 2) = 0:0. 

B(1,4) = -Krdot; 

B(2,2) co 1.0; 

B(2,4) = 0.0; 

B(3,2) = 0.0; 

B(3,4) = (M*Xg)-Yrdot; 
B(4,2) = 0.0; 

B(4,4) = Izz-Nrdot; 


time(1) = 0.0; 

p(1) = 0.0; pdot(1) = 0.0; 
phi(1)=1.0*pi/180; — phidot(1) = 0.0; 
psi(1) = 0.0; psidot(1) = 0.0; 
v(1) = 0.0; vdot(1) = 0.0; 
m1) = 0.0; rdot(1) = 0.0; 
x(1) = 0.0; y1)= 0.0; 
xdot(1) = 0.0; ydot(1) = 0.0; 
tfinal = 1; 


% THIS SECTION ALLOWS FOR RANDOM DISTURBANCES IN SWAY 
% AND YAW FOR THE LINEAR EQUATIONS AS WELL. MATLAB 

% ALLOWS FOR EITHER A UNIFORM OR NORMAL DISTRIBUTION 
% OF RANDOM NUMBER. 


rand(‘uniform’) 
vdotdist = rand(1,(tfinal/deltat)+1); 
vdotdist = 0.05*(vdotdist - mean(vdotdist)); 
rdotdist = rand(1,(tfinal/deltat)+1); 
rdotdist = 0.04*(rdotdist - mean(rdotdist)); 


% COMMENCE EULER INTEGRATION 
for j = 2:(tfinal/deltat) + 1; 


py) = p(j-1) + pdotG-1)*deltat; 
phi(j) = phi(j-1) + phidot(j-1)*deltat; 
psij)=  psi(j-1) + psidotG-1)*deltat; 
v= v(j-1) + vdotG-1)*deltat; 
beta(j) = -atan(v(j)/U0); 

re r(j-1) + rdotG-1)*deltat; 


x(j) = x(j-1) + xdot(j-1)*deltat; 

y(j) = y(j-1) + ydot(j-1)*deltat; 

pdot(j) = C(1,1)*pG) + C(1,2)*phig) + C(1,3)*vq) + C47); 
phidot(j)= C(2,1)*pG) + C(2,2)*phiG) + C(2,3)*vG) + C(2,4)*rG); 
vdot(j) = C(3,1)*pQ) + C(3,2)*phiG) + C(3,3)*vG) + C(3,4)*rQ); 
rdot(j) = C(4,1)*pG) + C(4,2)*phiG) + C(4,3)*vG) + C(4,4)*rq); 
psidot(j) = r(j)*cos(phi(j)); 

xdot(j) = U0*cos(psi(j)) - vG)*sin(psiG))*cos(phiGj)); 

ydot(j) = U0*sin(psi(j)) + vG)*cos(psiGj))*cos(phi(j)); 

time(j) = (j*deltat) - deltat; 

end 

return 
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% THIS IS A NUMERICAL INTEGRATION PROGRAM TO CALCULATE 
% THE DRAG FORCES IN THE HORIZONTAL PLANE UTILIZING THE 
% TRAPEZOIDAL RULE. THE CALL TO ‘crossflowintmod’ IS 

% IDENTICAL, ONLY USING DIFFERENT VARIABLE NAMES AND 

% VALUES FOR THE MODIFICATION PORTION OF HAMMING’S 

% METHOD. 


function[SWAY,YAW]=crossflowint(tmp1,tmp2,X,HH) 


RHO = 1.94; CDy=0.35; SWAY=0.0; YAW =0.0; 
v=tmpl;r=tmp2; 


for k = 1:18; 
if abs(v+X(k)*r) < le-6 
UCF(k) = 0.0; 
else 
UCF(k) = (v+X(k)*r)abs(v+X(k)*r)); 
end 


CFLOW(k) = CDy*HH(k)*((v+X(k)*r)*2); 
SWAY 1(k) = CFLOW(k)*UCF(k); 
YAWI1(k) = CFLOW(k)*X(k)*UCF(k); 
end 


m= 1:17; 
SWAY = -0.25*RHO*sum((SWAY 1Gj)+SWAY1(jj+1)).*(XGj+1)-XGj))); 
YAW = -0.25*RHO*sum((YAW1(jj)+YAW1(Gj+1)).*(XGj+1)-XGj))); 
return 


EIGENVALUE CALCULATION PROGRAM; 
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% THIS PROGRAM DETERMINES THE FOUR ROOTS FOR THE ROLL 
% COUPLED EQUATIONS OF MOTION OVER A RANGE OF Xg. 


function[all_roots,xg] = findroots(U0,Xgmin,stepXg,Xgmax) 


W = 12000.0; Ixx=1760.0; Iyy=9450.0; Izz=10700.0; 
L=17.425; RHO=1.94; G=32.2; M = W/G; 
BUOY = W;; Zb=0.0; 

ZZ = 0.5*RHO*L‘2; 


% SWAY HYDRODYNAMIC COEFFICIENTS 
Ypdot = 1.27E-04*ZZ*L*2;  Yrdot = 1.24E-03*ZZ*L/2; 
Yvdot = -5.55E-02*ZZ*L;  Yp = 3.055E-03*ZZ*L; 
Yv = -9.31E-02*ZZ; 


CC = input(‘enter model configuration choice; either 1 for A or 2 for B’) 


if CG 

Yr = -3.500e-02*ZZ*L; 
else 

Yr = -5.940e-02*ZZ*L; 
end 


% ROLL HYDRODYNAMIC COEFFICIENTS 
Kpdot =-1.01E-03*ZZ*L‘*3; Krdot =-3.37E-05*ZZ*L‘3; 
Kvdot = 1.27E-04*ZZ*L*2; Kp =-1.10E-02*ZZ*L‘2; 
Kr = -8.41E-04*ZZ*L“2; Kv = 3.055E-03*ZZ*L; 


% YAW HYDRODYNAMIC COEFFICIENTS 
Npdot = -3.370E-05*ZZ*L‘3; Nrdot = -3.400E-03*ZZ*L‘3; 
Nvdot = 1.240E-03*ZZ*L‘2; Np =-8.405E-04*ZZ*L‘2; 
Nr =-1.640E-02*ZZ*L“%2; Nv = =-1.484E-02*ZZ*L; 


rows=(abs(Xgmax-Xgmin)/stepXg)+ 1; 
all_roots = zeros(rows,4); xg = zeros(rows,1); 


Zg = input(‘enter value for Zg’') 


a = (Ixx-Kpdot); b = (Kp*U0); c = (Zg*W)-(Zb* BUOY); 
d = (M*Zg)+Kvdot; e = Kv*U0; f = Krdot; 

g = (M*Zg*U0) + (Kr*U0); h = Ypdot + (M*Zg); 

P= Yor U0: j=M- Yvdot; k = Yv*U0; 

o = Npdot; u = Izz - Nrdot; p = Np*U0; 

m = (Yr*U0) - (M*U0); x= Ny U0; 


Xg = Xgmin:stepXg:Xgmax; 
Xb = input(‘enter either 0.0 or Xg for value of Xb); 


1 = (M*Xg) - Yrdot; q = (Xb.*BUOY)-(Xg*W); 
r = (M*Xg)-Nvdot; w = (Nr*U0)-(M*Xg*U0); 


A = (((a).*(4*u-l.*r)+(d).*(u*h+o.*1)-(f).*(r.*h-+0*)))); 
B = ((e).*(u*h+o.*1)-(d).*(u*i-w.*h+o.*m-p.*1)-(a).*... 


G.*w+k*u-].*x-r.*m)-(b).*G*u-l.*r)+(f).*(r.*1-x... 
.*h+o*k-p*j)-(g).*(r.*h+0*))); 
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C = ((a).*(k.* w-m.*x)+(b).*G.*w+k*u-].*x-r.*m)}+(c).*... 
(*u-l.*r)-(d).*(q.*]-w.*i+m.*p)-(e).*(u*i-w.*h+... 
o.*m-p.*])+(f).*(q.*j-x*it+p*k)+(g).*(r.*i-x.*h+o*k... 
aBeed)): 

D = ((g).*(q.49-x.41+p*k)-(f.*q*k}+(d.*q.*m)}-(e).*(q.*L... 
-w.*i+m.*p}-(c).*G.*w+k*u-l.*x-r.*m)-(b).*(k.*w-m. *x)); 

FE = ((c).*(k.*w-m.*x)+(e.*g.*m)-(g.*q.*k)); 

Z=— (Aes CD EB}; 
j= l:rows; poly(j,=z(j,:); rts = ‘roots(poly(j,:))’; 
for jj) = 1:rows; 


all_roots(jj,1:4) = eval(rts)’;; 
end 


xg=Xg ; 


end; 


FIRST AND SECOND ORDER PERTURBATION 
APPROXIMATION COMPUTER PROGRAM: 
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% THIS PROGRAM DEVELOPS THE 1ST & 2ND ORDER PERTURBATION 
% APPROXIMATION SOLUTIONS FROM SECTION IV 


function[xlocation,perturbroot1,perturbroot2,xg] =perturbanalysis(U0) 


W = 12000.0; Ixx = 1760.0; Iyy = 9450.0; Izz=10700.0; 
L= 17.425; RHO=1.94; G=32.2; M = W/G; 
BUOY = W;; Zb=0.0; Xb = 0.0; 


ZZ = 0.5*RHO*L“2; 


% SWAY HYDRODYNAMIC COEFFICIENTS 
Ypdot = 0.0; Yrdot= 1.24§-03*ZZ*L“2; 
Yvdot = -5.55E-02*ZZ*L; Yp = 0.46, 
Yv =-9.31E-02*ZZ; 


CC = input(‘enter model configuration choice; either 1 for A or 2 for B’) 


i CC<2 
Yr = -3.500e-02*ZZ*L; 


else 
Yr = -5.940e-02*ZZ*L; 
end 


% ROLL HYDRODYNAMIC COEFFICIENTS 


Kpdot = -1.01E-03*ZZ*L‘3; Krdot = 0.0; 
Kvdot = 0.0; Kp = =-1.10E-02*ZZ*L‘2, 
Kr = -8.41E-04*ZZ*L“2; Kv = 3.055E-03*ZZ*L; 

% YAW HYDRODYNAMIC COEFFICIENTS 
Npdot = 0.0; Nrdot =-3.400E-03*ZZ*L‘43; 
Nvdot = 1.240E-03*ZZ*L‘2; Np = 0.0; 
Nr  =-1.640E-02*ZZ*L“2; Nv = =-7.42E-03*ZZ*L; 


Zg = input(‘enter value for Zg’) 


% THIS PERTURBATION PROGRAM SOLVES FOR THE SOLUTION 
% ABOUT Xg=0 


Xg = 0.0; 
index = 1; 


for Xg = -1.0:0.05:0.2; 


a = (Ixx-Kpdot); b = (Kp*U0); c = (Zg*W)-(Zb* BUOY); 
d = (M*Zg); e = Kv*U0; f = Krdot; 

g = (M*Zg*U0)+(Kr*U0); h= Ypdot+(M*Zg); i= Yp*U0; 

j= M-  Yvdot; k = Yv*U0; ] = (M*Xg) - Yrdot; 
m= (Yr*U0Q) - (M*U0); o = Npdot; p= Np*U0; 


q = (Xb*BUOY)-(Xg*W); = r = (M*Xg)-Nvdot; u = Izz - Nrdot; 
w =(Nr*U0)-(M*Xg*U0); x = Nv*U0; 


Al = G*u)-(/*r); Bl = -((u*k)+G*w)-(m*r)-(1*x));  C] = (k*w)-(x*m); 
Ar=a; Br=-b; Cre=c; 


K1 = d*u*e + d*Nvdot*Kr*U0; 
K2 = d*((Nr*#Kv*U0%2)-(Nv*Kr*U02)); 


alpha =-M*M*Zg*Kr*U0; 
beta = -d*((W*Yrdot)+(M*Kv*U0*U0)); 

gamma = (d*W*U0*(Yvdot-Yr))+(W*U0*(Kr* Yvdot-Kr*M-Kv*Yrdot)), 
delta =g*Yv*W*U0 - m*Kv*W*U0; 


%PQRS &T SOLVE THE 4TH ORDER POLYNOMIAL FOR a0 


P = Ar*Al; 

Q = Ar*Bl + Br*Al +K1; 

R = Ar*Ci +Cr*Al +Br*B! +K2; 
S = Cr*Bl + Cl*Br; 

LSet BE 


pa 1? 

B = Q+alpha*Xg; 
C = R+beta*Xg; 

D = S+gamma*Xg; 
E = T+delta*Xg; 


z=(ABCDE]; temp 1(index, 1:4) = roots(z)’; 
xe(index,1) = Xg; index = index + 1; end 
Za Go Ro Ti: temp2(1,1:4) = roots(zz)’; 


a0 = temp2(1,4) 
% THIS PART SOLVES FOR al IN THE 1ST ORDER PERTURBATION 


num1 = -(alpha*(a0‘%3)+beta*(a0%2)+gamma*a0+delta); 
denl = (4*P*(a0%3)+3*Q*(a0%2)4+2*R*(a0}S); 
al = num1i/denl 


% THIS PART SOLVES FOR a2 IN THE 2NB ORDER PERTURBATION 


num2a= -(a1%2*(6*P*a0*a0+Q*a0+R)+al*(gamma+2*beta*a0+... 
3*alpha*a0*a0)); 
a2 = num2a/den1 


index = 1; 
for X = -0.8:0.04:0.4; 


% PERTURBROOT1 IS 1ST ORDER PERTURBATION APPROXIMATION 
perturbroot1(index) = a0 + al*X; 


% PERTURBROOT2 IS 2ND ORDER PERTURBATION APPROXIMATION 
perturbroot2(index) = a0 + al*X + a2*X*X; 


xlocation(index) = X; 
index = index+1; 
end 

return 
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